Persistence of Tembusu Virus in Culex tritaeniorhynchus in Yunnan Province, China

The Tembusu virus (TMUV), a member of the Flaviviridae family, can be transmitted via mosquitoes and cause poultry disease. In 2020, a strain of TMUV (YN2020-20) was isolated from mosquito samples collected in Yunnan province, China. In vitro experiments showed that TMUV-YN2020-20 produced a significant cytopathic effect (CPE) in BHK, DF-1, and VERO cells, while the CPE in C6/36 cells was not significant. Phylogenetic analysis revealed that the strain belonged to Cluster 3.2 and was closely related to the Yunnan mosquito-derived isolates obtained in 2012 and the Shandong avian-derived isolate obtained in 2014. Notably, TMUV-YN2020-20 developed five novel mutations (E-V358I, NS1-Y/F/I113L, NS4A-T/A89V, NS4B-D/E/N/C22S, and NS5-E638G) at loci that were relatively conserved previously. The results of this study demonstrate the continuous circulation and unique evolution of TMUV in mosquitoes in Yunnan province and suggest that appropriate surveillance should be taken.


Introduction
The Tembusu virus (TMUV), which is a single, positive-stranded RNA virus, is a member of the Flaviviridae family. In Malaysia, Culex tritaeniorhynchus was reported as the first source of TMUV in 1955 [1]. Since then, other Culex species, including Culex gelidus and Culex vishnui, have also been found to carry the virus [2,3]. In 2000, researchers found TMUV could infect vertebrates and cause avian disease. TMUV was isolated from samples of sick chickens from a broiler flock in Malaysia and characterized by symptoms such as encephalitis and growth retardation [4].
In 2010, TMUV was found to spread among ducks in southeast China and caused a substantial economic effect [5]. The farming sector was seriously threatened when similar epidemics were discovered in duck populations in Southeast Asia, notably Malaysia and Thailand [6,7]. With a morbidity rate of 90% to 100% and a mortality rate between 5% and 30% in infected duck flocks, TMUV infection can have clinical indications in ducks, such as growth retardation, loss of appetite, paralysis, and decreased egg production [5,8,9]. TMUV not only affects poultry (such as chickens, ducks, and geese), but also infects birds (such as pigeons and sparrows) [10][11][12] with similar clinical symptoms.
Mosquito-borne flaviviruses, such as Japanese encephalitis virus (JEV), dengue virus (DENV), and Zika virus (ZIKV), can infect humans and cause large-scale disease outbreaks [13,14]; nevertheless, no cases of TMUV infection in humans or other mammals have been reported. However, some studies have shown that intracerebral inoculation of BALB/c mice with TMUV causes severe neurological symptoms and even death [15]. Furthermore, epidemiological studies have shown that IgG antibodies or neutralizing antibodies to TMUV have been detected in the sera of some healthy people in China and Thailand [16,17]. These findings imply that TMUV can cause considerable economic losses in the poultry farming industry and pose a potential threat to public safety.
The TMUV genome contains 11,000 nucleotides that encode 10 proteins, including 3 structural proteins: capsid (C), pre-membrane (prM), and envelope (E); and 7 nonstructural (NS) proteins: NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5 [18]. The E protein is frequently used for phylogenetic analysis of genetic variation and molecular epidemiology as it is the main surface protein of the viral particle. It also plays a crucial role in mediating virus receptor binding and virus cell membrane fusion and frequently undergoes genetic mutations under the evolutionary pressure of the host immune system [19]. TMUV strains are currently divided into 3 clusters (Cluster 1, 2, and 3) based on phylogenetic analysis of the open reading frame (ORF) region and the E gene [20,21]. Cluster 2 is currently the dominant cluster in China and Cluster 3 encompasses the majority of mosquito-derived isolates [22].
Yunnan province is located at the southwest border of China, where a considerable number of arboviruses have been isolated [23]. To monitor arboviruses in Yunnan province, we collected mosquitoes from July to August 2020 in Yunnan province and successfully isolated a strain of TMUV from Culex tritaeniorhynchus. The new strain belongs to Cluster 3.2 and has five novel amino acid mutations. Our findings reveal that TMUV is still spreading in mosquitoes in Yunnan province.

Mosquito Collection
The mosquito samples were collected from July to August 2020 in Nansa village, Yuanjiang county, Yuxi city, Yunnan province using UV mosquito traps (Kung Fu Xiaoshuai, 12 V, 300 mA; Wuhan Jixing Medical Technology Co., Ltd., Hubei, China) in pig and cattle pens. The traps were set between 18:00 p.m. and 6:00 a.m. The collected mosquito samples were morphologically characterized and mosquitoes of the same species were maintained in liquid nitrogen in pools of 50-100 until they were tested in the laboratory [24].

Virus Isolation and RT-PCR Analysis
The collected mosquito samples were ground as previously reported [25]. The ground supernatants were inoculated into culture plates containing monolayers of BHK-21 and C6/36 cells grown at 37 • C, 5% CO 2 , and 28 • C, respectively. Cell supernatants were collected and inoculated into newly cultured BHK-21 and C6/36 cells for a third passaging after 1 observation cycle (6-7 days). The positive isolates that could produce cytopathic effects (CPEs) were collected after three generations for identification.
The total RNA was extracted from virus-infected cells according to the instructions of the QIAamp Viral RNA Mini Kit (Qiagen, Valencia, CA, USA). Reverse transcription was performed with Ready-to-Go TM You-Prime First-Strand Beads (GE HealthCare, Little Chalfont, Buckinghamshire, UK) and random primers. Flavivirus genus-specific primers (F: 5 -TACCACATGATGGGAAAGAGAGAGAGAA-3 , R: 5 -GTGTCCCAGCCGGCGGTGTC ATCAGC-3 ) [26] were used for the RT-PCR detection, and the reaction program was set up in accordance with the instructions on the PrimeScriptTM One Step RT-PCR Kit Ver. 2 (Dye Plus) (TaKaRa Bio Inc, Shiga, Japan) reagent package. The reaction parameters were as follows: 94 • C for 5 min; denaturation at 94 • C for 30 s, annealing at 55 • C for 30 s, extension at 72 • C for 30 s for 35 cycles; final extension at 72 • C for 10 min. The intended target product size was 310 bp and the RT-PCR amplification products were sequenced using Sanger sequencing.

Viral Titer Determination
The virus isolate was propagated in a BHK-21 cell and a standard plaque assay was used to calculate the virus titers. Briefly, a 6-well culture plate was filled with a succession of 10-fold dilutions of the virus. Each well received 4 mL of 1.1% methylcellulose-MEM semisolid media with 2% FBS after incubation. After 4 days of incubation, we noticed the plaques clearly developing. The plaque formation units and virus titers were calculated after staining cells with crystalline violet solution.

Virus Growth Kinetics
As previously described, BHK-21, DF-1, VERO, and C6/36 cells were infected with the virus at a multiplicity of infection (MOI) of 0.1, respectively. After 1 h of incubation, the cells were gently washed with PBS (Invitrogen, Carlsbad, CA, USA), then the cell culture medium was added, and the cells were placed under appropriate conditions. The cell and supernatant samples were collected at 0, 12, 24, 36, 48, 60, 72, and 96 hpi and stored at −20 • C. A standard plaque assay was used to calculate viral titers.

Virus Genome Sequencing
For high-throughput sequencing, RNA from the flavivirus positive isolate's third generation cell culture supernatant was used. Sequencing libraries were created using the Trio RNA-SeqTM library preparation kit guidelines (Tecan Genomics, Männedorf, Switzerland). The libraries were filtered with AMPure XP Beads (Beckman Coulter, Brea, CA, USA) and assessed for library quality with the Qubit4.0 Fluorometer (Life T Technologies, Grand Island, NY, USA). The library was then sequenced on an Illumina Novaseq 6000 platform, resulting in 150 bp paired-end reads. Offline data were filtered based on Q30 values and read lengths for further data analysis. The sequencing data were spliced using the reference genome NC015843 sequence as a template after low-quality reads were eliminated using CLC Genomics Workbench. Moreover, we sequenced the original mosquito pool using the Sanger method. Primers were acquired from previous articles (Supplement Table S1) [27]. The sequence was submitted to NCBI (accession number: OQ238827).

Sequence and Phylogenetic Analysis
TMUV sequences with temporal, host, and geographic information were acquired from GenBank, and the newly sequenced sequences from this work were added to them, resulting in a dataset of 148 TMUV sequences for analysis (Supplement Table S1). Multiple sequence alignment analysis was performed using MAFFT software (v7.450) [28]. MegAlign (DNAStar, Madison, WI, USA) was used to analyze viral sequences for nucleotide and amino acid homology [29]. A maximum likelihood (ML) tree reconstruction of 148 TMUV genome sequences was performed using IQ-TREE software (v1.6.12), and the best-fit alternative model for the TMUV genome determined by ModelFinder software was GTR + F + R3 [30]. The final tree files were visualized and annotated using the Chiplot online website (www.chiplot.online, accessed on 15 November 2022). Python scripts were used for mutation analysis in the gene sequences of the TMUV.

Bayesian Phylogenetic Reconstruction and Time to the Most Recent Common Ancestry (tMRCA)
The maximum clade credibility (MCC) tree was created using Beast software (v1.10.4) based on 148 full-length E-gene sequences. The alternative GTR + G + I model was employed in this investigation, with the strict clock was selected as the clock type and the Bayesian skyline as the tree prior, so Markov chain Monte Carlo (MCMC) was 200,000,000 and data gathering was conducted every 20,000 generations. To ensure that the estimated sampling size (ESS) were greater than 200, the results were analyzed visually using Tracer software (v1.7.1) [31]. The TreeAnnotator tool in the Beast package was used to process the data. The tree with the highest posterior probability (PP) with a 10% burn-in was chosen to represent a target tree. FigTree software 1.4.4 was used to visualize and comment on the final tree.

Mosquito Collection
A total of 2723 Culex tritaeniorhynchus and 2273 Anopheles sinensis were collected from Yunnan province, China, between July and August 2020.

Virus Isolation and Identification
The mosquitoes were combined into pools of 50-100. A total of 94 mosquito pools were homogenized and inoculated in parallel with the BHK-21 and C6/36 cells. One pool from the BHK-21 cells and six pools from the C6/36 cells showed CPE.
RT-PCR was used to analyze positive isolates using the Flavivirus genus primers and the results showed that only the pool from the BHK-21 cells was positive. This positive isolation derived from Culex tritaeniorhynchus was identified as TMUV by NCBI-BLAST for the Sanger sequencing of RT-PCR production and was named TMUV-YN2020-20.

Virus Growth Kinetics
Mosquito (C6/36), avian (DF-1), and mammalian (BHK-21, VERO) cell lines were used to compare the susceptibility to TMUV-YN2020-20 ( Figure 1a). The findings revealed that BHK-21 and DF-1 cells exhibited CPE at 48 hpi, and the cells shrank and gradually shed at 96 hpi. The VERO cell showed cell rounding at 96 hpi, whereas the C6/36 cells did not display significant CPE at 96 h. All four infected cell lines showed replication and proliferation ( Figure 1b). The replication of this strain in the BHK-21 and DF-1 cells showed an increasing trend, then decreased within 96 h. Viral titers peaked at 1 × 10 7.88 PFU/mL in DF-1 cells and at 1 × 10 7.49 PFU/mL in BHK cells at 60 hpi, in the later stages of infection due to cell death, the viral titer decreased until 96 hpi. In the VERO and C6/36 cells, the viral replication increased with the increasing infection time up to 96 h. At 96 hpi, the viral titers were 1 × 10 5.98 PFU/mL and 1 × 10 6.52 PFU/mL in the VERO and C6/36 cells, respectively.

Phylogenetic Analysis
The genome sequences of TMUV-YN2020-20 were obtained by high-throughput sequencing, and the total genome length for this strain was 10,994 nt, containing the 5′-UTR (94 nt), ORF region (10,275 nt/3425 aa), and 3′-UTR (625 nt). The result of Sanger sequencing was consistent with the high-throughput sequencing.

Phylogenetic Analysis
The genome sequences of TMUV-YN2020-20 were obtained by high-throughput sequencing, and the total genome length for this strain was 10,994 nt, containing the 5 -UTR (94 nt), ORF region (10,275 nt/3425 aa), and 3 -UTR (625 nt). The result of Sanger sequencing was consistent with the high-throughput sequencing.
Additionally, we identified 66 cluster-specific amino acids between Cluster 1, Cluster 2, and Cluster 3 TMUVs, most of which were located in the E protein, followed by NS1, NS2A, and NS5 proteins ( Table 2). By analyzing the amino acid variability in ORF

Discussion
A TMUV strain, TMUV-YN2020-20, was isolated from Culex tritaeniorhynchus derived from cattle farms in Yunnan province, China. This strain was replicated in the BHK-21, DF-1, VERO, and C6/36 cells, and had significant CPE in all cells except C6/36. The analysis of the viral genome revealed that it belonged to Cluster 3.2, which was most closely related to the YN12193 isolate at the nucleotide/amino acid level, as well as at the genetic evolutionary level. In addition, this strain evolved five novel amino acid mutations.
TMUV-YN2020-20 does not produce significant CPE in C6/36 cells, and this phenomenon also existed in the study of TMUV-TP1906 strain [32]. However, these two strains could replicate in C6/36 cells. A reasonable explanation is that mosquito-derived TMUV has evolved over a long period to adapt to mosquitoes and mosquito-derived cell lines, and then is transmitted by mosquitoes. Culex mosquitoes, such as Cx. pipiens and Cx. tritaeniorhynchus, are the main vectors of TMUV in nature [27,33]. It has not been determined whether Aedes mosquitoes can act as vectors of TMUV; however, in studies on the efficacy of mosquito transmission, it was found that Aedes albopictus can be infected with TMUV [34]. Additionally, A. albopictus can act as an effective vector of virus transmission when the titer of the infecting virus reaches a certain level, indicating the importance of mosquitoes in TMUV transmission.
Phylogenetic and population dynamic analyses pointed out that Cluster 1 contained avian-derived TMUV strains from Thailand and Malaysia; most of the avian-derived TMUV strains from China were located in Cluster 2, while Cluster 3 contained mosquito-derived TMUV strains. Furthermore, these mosquito-derived strains cluster with the avian-derived strains, indicating that TMUV has the potential for mosquito-to-avian transmission. The Yunnan mosquito-derived TMUVs (YN12115, YN12193, and YN2020) belong to Cluster 3.2 and these Yunnan strains share the most recent common ancestor, suggesting that they have the same evolutionary origin. TMUV is currently thought to survive only in mosquitoes or avian species; however, no cases of the virus infection in diseased ducks were found in Yunnan province until 2019 [35]. Additionally, TMUV was last found in mosquitoes from Yunnan in 2012. We consider that mosquitoes in Yunnan province continue to carry the virus due to the similarity between the mosquito-derived TMUVs found in 2012 and 2020. Furthermore, a conjecture is that TMUV may only circulate within mosquitoes before duck infection, or in a more covert manner in mosquito-duck; nonetheless, either scenario requires mosquitoes to play a key role. In addition, we could not analyze the more precise association between mosquito-derived and duck-derived TMUVs in Yunnan province since the sequences of duck-derived TMUVs are not yet available in the database. However, TMUV isolated from the Shandong avian (duck) was present in Cluster 3.2; therefore, it is likely that these mosquito-derived viruses have the ability to transmit to ducks. Long-term continuous monitoring is necessary to better investigate the circulation of TMUV in natural sources. In brief, the re-isolation of TMUV from mosquitoes in Yunnan province after an interval of 8 years demonstrates the persistence of the virus in this area.
There are two main methods of mosquito-borne virus transmission: one is horizontal transmission, including mosquito-animal (mosquito bite transmission) and mosquitomosquito (sexual transmission), the other is vertical transmission (ovarian transmission) [36]. Mosquito-animal is the general mode of transmission, and several relevant arboviruses, such as JEV, ZIKV, and DENV, use this mode [37,38]. Mosquito-transmitted TMUV is also considered to rely on this mode [39]. Sexual transmission is a specific type of transmission by which female mosquitoes transmit the virus to male mosquitoes that do not suck blood to expand the ecological range of the virus; ZIKV and others have been found to use this method [40]. Vertical transmission is a mechanism by which viruses can maintain circulation in nature under unfavorable conditions [41]. Members of the family Flavivirus, such as DENV [42] and ZIKV [36], have been found to be transmitted to offspring via eggs in Aedes aegypti and Aedes albopictus. No study has confirmed that TMUV can be sexually transmitted. In laboratory studies, it was concluded that duck-derived TMUV DK/TH/CU-1 cannot be transmitted vertically in mosquitoes [43]. Whether mosquito-derived TMUV can be transmitted vertically in mosquitoes and its mode of transmission in birds require additional studies.
TMUV-YN2020-20 has five new amino acid mutations (V358I in the E protein, Y/F/I113L in the NS1 protein, T/A89V in the NS4A protein, D/E/N/C22S in NS4B protein, and E638G in the NS5 protein). Amino acid mutation may contribute to the adaptive evolution of the virus to vector [44]. It has been shown that the glycosylation site at amino acid 154 of most flavivirus E proteins can facilitate virus transmission by overcoming the tissue barrier or antiviral response in mosquitoes [45]. Furthermore, amino acid residues at 304 and 367 of the E protein are associated with the virulence and pathogenicity of the virus [46,47]. In addition to the E protein, mutations at amino acid residues 543 and 711 of the NS5 protein affect viral replication and infectivity in the host [48]. The effects of these novel mutations on the function of viral proteins, viral replication, and virulence of the virus require further investigation.
In conclusion, a newly discovered TMUV isolate was analyzed for cellular adaptation and gene sequence characteristics, emphasizing the diversity of TMUV genetic evolution. At this stage, most studies on the pathogenicity of TMUV have focused on Cluster 2, while the remaining genotypes have been less studied, and mosquito-to-avian transmission patterns have not been studied in depth. It is necessary to conduct continuous monitoring of mosquitoes to further investigate the epidemiology of TMUV in Yunnan, China.